
library(sna)
library(lme4)

data1<-read.csv("Outdset.csv")
dataNOIN<-subset(data1, data1$state!="IN")
data1$committeeorigin<-ifelse(is.na(data1$committeeorigin)==TRUE, "No Committee", data1$committeeorigin)


#.75 sds
#data1$wksupport<-data1$wksupport75
#data1$stsupport<-data1$stsupport75
#data1$alteroutwk<-data1$alterinwk75
#data1$alteroutst<-data1$alterinst75

#1.25 sds
#data1$wksupport<-data1$wksupport125
#data1$stsupport<-data1$stsupport125
#data1$alteroutwk<-data1$alterinwk125
#data1$alteroutst<-data1$alterinst125



#Normalize Across Chambers for convergence#


data1$alterwk.center<-(data1$alteroutwk-mean(data1$alteroutwk))/sd(data1$alteroutwk)
data1$alterst.center<-(data1$alteroutst-mean(data1$alteroutst))/sd(data1$alteroutst)
data1$wk.center<-(data1$wksupport-mean(data1$wksupport))/sd(data1$wksupport)
data1$st.center<-(data1$stsupport-mean(data1$stsupport))/sd(data1$stsupport)

model.2<-lmer(outofcomm~insadv+tenure+cosponsor+party+alterwk.center+alterst.center+wk.center+st.center+(1|state), family=binomial (link=logit), data=data1)
model.2

model.2alter<-lmer(outofcomm~insadv+tenure+cosponsor+party+alterwk.center+alterst.center+(1|state), family=binomial (link=logit), data=data1)
model.2alter

model.3<-lmer(outofcomm~insadv+tenure+cosponsor+party+wk.center+st.center+(1|state), family=binomial (link=logit), data=data1)
model.3

model.4<-lmer(outofcomm~insadv+tenure+cosponsor+party+wk.center+alterwk.center+(wk.center:alterwk.center)+(1|state), family=binomial (link=logit), data=data1)
model.4  

model.5<-lmer(outofcomm~insadv+tenure+cosponsor+party+st.center+alterst.center+(st.center:alterst.center)+wk.center+alterwk.center+(wk.center:alterwk.center)+(1|state), family=binomial (link=logit), data=data1)
model.5

model.g<-glm(outofcomm~insadv+tenure+cosponsor+party+alterwk.center+alterst.center+wk.center+st.center+as.factor(state), family=binomial (link=logit), data=data1)
vif(model.g)   

###########DO YOUR SELECTION MODEL IN STATA#############


#Predicted Probabilities

#Non Conditional Model
library(mvtnorm)
sims<-10000
pp.sim<-rnorm(sims, fixef(model.2)[8], vcov(model.2)[8,8])
calc<-50
z<-seq(-3.6, 5.6, length=calc)
lo<-numeric(calc)
hi<-numeric(calc)
pe<-numeric(calc)
for(i in 1:calc){
pp<-numeric(sims)
for(j in 1:sims){
pp[j]<-pp.sim[j]*z[i]
}
pe[i]<-mean(invlogit(pp))
lo[i]<-quantile(invlogit(pp), 0.001)
hi[i]<-quantile(invlogit(pp), 0.999)
}

plot(z, pe, type="l", ylim=c(0.3,0.7), xlab="Secondary Connections from Weak Ties", ylab="Probability of Bill Survival", main="Bill Success in Committee and Sponsor Secondary Ties")
lines(z, lo, type="l", lty=2, ylim=c(0,1)) 
lines(z, hi, type="l", lty=2, ylim=c(0,1)) 
abline(0.5, 0)
legend(-2, 0.7, c("Point Estimates", "95% CI's"), col=c("black", "black"), lty=c(1, 2))

#Combined Model

#Marginal Effects Graph                                    

pdf("MargFXDirectWeak.pdf")
L<-seq(-1.1, 2.6, by=0.01)
margfx3<-1.64*(sqrt(vcov(model.5)[8,8]+vcov(model.5)[11,11]*L^2+2*L*vcov(model.5)[8,11]))
plot(L, fixef(model.5)[8]+fixef(model.5)[11]*L,  type="l", col="black", ylim=c(-0.25, 0.6), lwd=1, cex.lab=1.25, cex.axis=1.25,cex.main=1.25, ylab="Marginal Effect of Weak Ties", xlab="Number of Secondary Ties (Standardized)", main="Effect of Weak Ties as Secondary Ties Vary")
lines(L, fixef(model.5)[8]+fixef(model.5)[11]*L+margfx3, type="l", lty=2, col="black", lwd=1)
lines(L, fixef(model.5)[8]+fixef(model.5)[11]*L-margfx3, type="l", lty=2, col="black", lwd=1)
abline(0,0, lwd=1)
dev.off()

pdf("MargFXDirectStrong.pdf")
L<-seq(-1.2, 4.1, by=0.01)
margfx3<-1.64*(sqrt(vcov(model.5)[6,6]+vcov(model.5)[10,10]*L^2+2*L*vcov(model.5)[6,10]))
plot(L, fixef(model.5)[6]+fixef(model.5)[10]*L,  type="l", col="black", ylim=c(-0.25, 0.6), lwd=1, cex.lab=1.25, cex.axis=1.25,cex.main=1.25, ylab="Marginal Effect of Stong Ties", xlab="Number of Secondary Ties (Standardized)", main="Effect of Strong Ties as Secondary Ties Vary")
lines(L, fixef(model.5)[6]+fixef(model.5)[10]*L+margfx3, type="l", lty=2, col="black", lwd=1)
lines(L, fixef(model.5)[6]+fixef(model.5)[10]*L-margfx3, type="l", lty=2, col="black", lwd=1)
abline(0,0)
dev.off()






#3d?? Need a matrix for dep var
library(arm)
q<-seq(-3.9, 3.5, length=10)
r<-seq(-3.6, 5.6, length=10)
model<-function(a,b){
invlogit(fixef(model.5)[8]*a+fixef(model.5)[9]*b+fixef(model.5)[11]*a+fixef(model.5)[11]*b)
}

z<-outer(q,r,model)

e<-seq(-2.2, 4.8, length=10)
f<-seq(-2.5, 4.2, length=10)
model<-function(a,b){
invlogit(fixef(model.5)[6]*a+fixef(model.5)[7]*b+fixef(model.5)[10]*a+fixef(model.5)[10]*b)
}

z1<-outer(e,f,model)

#pdf("3dplot1.pdf")
persp(q,r,z, ticktype="detailed", shade=0.2, theta=140, phi=0, expand=0.8, lty=1, lwd=1, col="white", cex.lab=1.25, cex.axis=1.2, cex.main=1.2, xlab="Secondary Ties", ylab="Direct Weak Ties", zlab="Pr(Bill Survival in Committee)", zlim=c(0.2, 0.85), main="Probability of Survival as Weak Direct and Secondary Ties Change")
#legend (x=-0.38, y=0.35, legend=c("Weak Ties", "Strong Ties"), fill=TRUE, col="00000088")
#dev.off()


#pdf("3dplot2.pdf")
persp(e,f,z1, ticktype="detailed", shade=0.2, theta=140, phi=0, col="#00000088", lty=1, lwd=1, expand=0.8, cex.lab=1.25, cex.axis=1.2, cex.main=1.2, xlab="Secondary Ties", ylab="Direct Strong Ties", zlab="Pr(Bill Survival in Committee)", zlim=c(0.2, 0.85), main="Probability of Survival as Strong Direct and Secondary Ties Change")
#legend (x=-0.38, y=0.35, legend=c("Weak Ties", "Strong Ties"), fill=TRUE, col="00000088")
#dev.off()
#00000088

###################################
#Difference in Marginal Effects Plots
library(mvtnorm)
L<-seq(-1.2, 2.6, by=0.01)

 cond.sim <- function(ind1, coeff, vcv){
                ind2 <- (1:length(coeff))[-ind1]
                sig11 <- vcv[ind1, ind1]
                sig12 <- vcv[ind1, ind2]
                sig21 <- vcv[ind2, ind1]
                sig22 <- vcv[ind2, ind2]
                mu1 <- coeff[ind1]
                mu2 <- coeff[ind2]
                mu.bar <- mu1
                sig.bar <- as.matrix(sig11 - (sig12%*%solve(sig22))%*%sig21)
return(list(mu.bar, sig.bar))
}

object.1 <- cond.sim(c(6, 10, 8, 11), fixef(model.5), vcov(model.5))


yy<-rmvnorm(10000, object.1[[1]], object.1[[2]])


ped<-(mean(yy[,3])+mean(yy[,4])*L)-(mean(yy[,1])+mean(yy[,2])*L)
ul<-(quantile(yy[,3], 0.95)+quantile(yy[,4], 0.95)*L)-(quantile(yy[,1], 0.95)+quantile(yy[,2], 0.95)*L)
ll<-(quantile(yy[,3], 0.05)+quantile(yy[,4], 0.05)*L)-(quantile(yy[,1], 0.05)+quantile(yy[,2], 0.05)*L)

pdf("DiffMFX.pdf")
plot(L, ped, type="l", ylim=c(0, 0.5), xlab="Range of Weak Ties (Standardized)", ylab="Difference between Marginal Effects", main="Difference between Weak and Strong Ties' Marginal Effects")
lines(L, ul, lty=2)
lines(L, ll, lty=2)
abline(0, 0, lwd=2)
legend(-1, 0.5, c("Difference in Marginal Effects","95% Confidence Intervals for Difference"), lty=c(1,2))
dev.off()

####Pred Prob Plots Conditional
library(arm)
LL<-seq(-1.2, 2.29, by=0.001)
object.2<-cond.sim(c(8,11), fixef(model.5), vcov(model.5))
yy2<-rmvnorm(10000, object.2[[1]], object.2[[2]])
plot(LL, x, type="l", ylim=c(0.4,0.65))
lines(LL, x2, col="red")
pp<-invlogit((mean(yy2[,1])+mean(yy2[,2])*-1.1)*LL)
pp2<-invlogit((mean(yy2[,1])+mean(yy2[,2])*2.6)*LL)
lo<-invlogit((quantile(yy2[,1], 0.05)+quantile(yy2[,2], 0.05)*-1.1)*LL)
lo2<-invlogit((quantile(yy2[,1], 0.05)+quantile(yy2[,2], 0.05)*2.6)*LL)
hi<-invlogit((quantile(yy2[,1], 0.95)+quantile(yy2[,2], 0.95)*-1.1)*LL)
hi2<-invlogit((quantile(yy2[,1], 0.95)+quantile(yy2[,2], 0.95)*2.6)*LL)

plot(LL, pp, ylim=c(0.4,0.65), type="l")
lines(LL, pp2, col="red")

#Contour plots?
library(arm)
q<-seq(-3.9, 3.5, length=10)
r<-seq(-3.6, 5.6, length=10)
model<-function(a,b){
invlogit(fixef(model.5)[8]*a+fixef(model.5)[9]*b+fixef(model.5)[11]*a+fixef(model.5)[11]*b)
}

z<-outer(q,r,model)

e<-seq(-2.2, 4.8, length=10)
f<-seq(-2.5, 4.2, length=10)
model<-function(a,b){
invlogit(fixef(model.5)[6]*a+fixef(model.5)[7]*b+fixef(model.5)[10]*a+fixef(model.5)[10]*b)
}

z1<-outer(e,f,model)

filled.contour(q,r,z, nlevels=60, col=grey(80:1/80), zlim=c(0.3,0.7), xlab="Weak Ties (Normalized)", ylab="Secondary Ties from Weak Ties", main="Predicted Probability of Bill Survival")

filled.contour(e,f,z1, nlevels=60, col=grey(80:1/80), zlim=c(0.3,0.7), xlab="Strong Ties (Normalized)", ylab="Secondary Ties from Strong Ties", main="Predicted Probability of Bill Survival")




